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Cycles  X  10 

Figure  2 

The  empiric  cumulative  and  the  distribution  F(a,8)  for  fatigue  life 
at  a  stress  of  26,000  psi. 
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Summary 


The  estimation  problem  Is  studied  for  a  new  two- 
parameter  family  of  life  length  distributions  which  has 
been  previously  derived  from  a  model  of  fatigue  crack 
growth.  Maximum  likelihood  estimates  of  both  parameters 
are  obtained  and  iterative  computing  procedures  are  given 
and  examined.  A  simple  estimate  of  the  median  life  is 
exhibited,  shown  to  be  consistent  and  then  compared, 
favorably,  with  the  maximum  likelihood  estimate.  More 
importantly  the  asymptotic  distribution  of  this  estimate 
is  shown  co  be  within  the  same  class  of  distributions 
as  the  observations  themselves.  This  model,  and  these 
estimation  procedures,  are  tried  by  fitting  this 
distribution  to  several  extensive  sets  of  fatigue  data 
and  then  some  comparisons  of  practical  significance 
are  made. 


1.  Introduction 


A  new  family  of  distributions  was  derived  in  [1]  from 
considerations  of  the  physical  behavior  of  fatigue  crack  growth 
under  cyclic  loading.  The  life  length  obtained  was  the  mathematical 
representation  of  the  number  of  cycles  needed  to  force  the  fatigue 
crack  to  exceed  a  critical  value. 

Let  us  denote  by  9  this  two-dimensional  parametric  family 
of  distributions  of  nonnegative  random  variables  defined  by 

F(t:ot,6)  -  9i[“  £(t/B)]  for  t  >  0  (1.1) 

where  a  >  0,  8  >  0  and 

£(t)  -  t 5  -  t"*5  (1.2) 

and  9?  is  the  distribution  function  of  the  standard  normal  variate. 

Of  course,  there  have  been  many  families  of  distributions  which 
have  been  suggested  as  candidates  for  the  theoretical  probability  law 
governing  fatigue  life.  Excellent  discussions  and  comparisons  of  many 
of  these  have  been  made  by  other  authors,  see  for  example  [3],  [5], 
and  [6]  and  the  references  given  there.  However,  the  main  intent  of 
this  study  is  to  investigate  the  parametric  estimation  problem  for 
the  family  defined  in  (1.1)  and  not  to  argue  its  particular  merits  in 
such  applications  over  other  families  of  distributions. 

In  this  paper,  we  derive  the  maximum  likelihood  estimates  of  the 
parameters  a  and  8  and  develop  some  iterative  numerical  procedures 
for  their  computation.  We  obtain  a  simple  estimate  of  8,  which  is 


shown  to  be  consistent  and  shown,  for  small  values  of  a ,  to  be 
virtually  the  same  as  the  maximum  likelihood  estimate.  In  all  cases, 
It  can  be  used  as  a  good  Initial  guess  for  the  iterative  procedures 
which  compute  the  maximum  likelihood  estimate.  We  also  obtain  the 
bias  and  variance  of  this  simplified  estimate  under  certain  limiting 
conditions.  Two  different  numerical  procedures  for  computing  the 
maximum  likelihood  estimate  of  g  are  proposed  and  their  behavior 
compared  by  applying  them  to  several  large  sets  of  observations  with 
different,  but  known,  values  of  the  parameters.  We  finally  present 
the  results  of  fitting  this  distribution  by  these  estimation 
procedures  to  some  rather  extensive  sets  of  fatigue  data  obtained 
from  metal  coupons  cycled  at  various  stress  levels. 

In  order  to  make  this  paper  self-contained  we  now  quote  without 
proof  some  results  given  in  [1],  which  will  be  used  subsequently. 

Theorem  1.1.  If  T  has  the  life  distribution  F(a,g)  in  & 
then  1/T  has  the  distribution  F(a,^-)  also  in  iF,  and  for  any 

P 

real  a  >  0  the  random  variable  aT  has  a  distribution  F(a,ag) 
in  Moreover, 


2 


ET  -  g(l  +  y~) 

(1.3) 

var(T)  -  (aS)2(l  + 

(1.4) 

and  if  Z  is  a  standard  normal  variate,  then 

2  -  /  2 

1  +  y  ZS  aZ/1  +  Z  (1.5) 

has  the  distribution  F(a, 1)  in 


2.  Maximum  Likelihood  Estimates  and  Their  Computation 


Now  we  can  state  the  primary 

Theorem  2.1.  If  t^,...,tn  is  a  sample  of  independent 
observations  each  distributed  by  T(a,P)  e  <7,  then  the  maximum 
likelihood  estimate  g  of  8  is  the  unique  positive  solution 
of  g(x)  ■  0  where  g  is  the  random  function  defined  by 

g(x)  -  x2  -  x[2r  +  K(x) ]  +  r [s  +  K(x)].  (2.1.3) 

Furthermore  r  <  6  <  s.  We  can  also  express  the  maximum  likelihood 
estimate  a  of  a,  defined  in  terms  of  6,  by 


(2.1.4) 


Proof.  Consider  the  density  for  t  >  0  obtained  from  (1.1) 


I 


Taking  the  natural  logarithm  of  the  joint  density  of  the  observations 

t.,...,t  we  obtain  the  likelihood 
1  n 

n  2 

L  s  -n  in  a  -  n  Hn  8  +  Y  {~h  «.n (2tt )  -  h  ct  £  (t  /g)  +  £n  ^  *  ( t . /g) } . 


ailL  2  l  £  2,  /on 

m  a.  -  —  2.  C  (t./g) 
n  3  a  n  i 


(2.2.1) 


aL  n  -2  £  i  £  t  C"(t,/g) 

3?  "  “  B  +  (“6)  ti^(ti/2H’(ti/B)  -  ^'(ti/6)  ‘ 


But  we  notice  that 


£(t)  «  /t  - 


£2(t)  -  t  +£ 


C*  (t) 


tVt-1* 


2£(t) 


(l-l/o 


(2.2.2) 


tCLU,  .  ,  +  1  tzl 

t'(t)  i  +  2t+l 


1  ±_ 

2  '  t+1  * 


Hence  by  substitution  we  have 


3L  -n  n  .  s_  n 

3B  "  2S  2a2g  d  ”  r'  +  K(B) 


(2.2.3) 


2^1  3L  =  2  , s  _  2a^. 

n  9B  r;  K(B)  ’ 


(2.2.4) 


By  equating  (2.2.1)  to  zero  we  then  have  the  equation 


2  s  ,  1  , 

a  =  ~  + - L 

6  r 


(2.3.1) 


and  by  equating  (2.2.4)  to  zero 
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a2  .  SL  .  £.  +  l&L 

6  r  K(g)  * 

Equating  (2.3.1)  and  (2.3.2)  and  simplifying  we  obtain 


(2.3.2) 


■fi  .  i  +  fig 

r  K(g)  * 


Substituting  from  (2.3.1)  for  a2  and  simplifying  we  have  that  the 
maximum  likelihood  estimation  of  g  is  the  solution  of  g(x)  -  0 
for  0  <  x  <  oo  (presuming  for  the  present  it  is  unique).  Now  we 
will  argue  that  g  is  unique  and  r  <  §  <  8.  Note  that  g(0)  - 
r[s+K(0)]  -  r(s+r)  >  0.  We  check  that  g(x)  +  -»  as  x  -►  ®,  by 
seeing  that  ^  -  1  and  [x-K(x)  ]  ,  -s  as  *  ,  ..  ^en 

■  x  -  K(x)  +  r  -  2r  + 

*  X  x 

so  g(x)/x  -►  -(s+r)  ae  x  -»■  <».  Now 

g'(x)  -  (x-r)U-K’(x)]  +  x  -  r  -  K(x)  (2. 


(2.4.1) 


and  one  sees  that 


K'(x)  ■  K2(x)  ^  Z(x+ti)-2, 


(2.4.2) 


and  we  know,  since  (Elx|V/v  is  nondecreasing  in  v  for  any  r.v.  X, 
we  have  K*  (x)  >  1.  Therefore,  since  x  -  K(x)  is  decreasing  we  see  that 
K(x)  >  x-r  for  x  >  0  and  we  have  g’(x)  <  0  for  x  >  r  with  at  most 
one  change  of  sign.  Thus  we  have  shown  g  unique  now  calculate 


g(r)  -  r(s-r) , 


g(s)  -  (s-r) [s-K(s) ] . 


t 
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Again  by  the  argument  above  we  know  that  s  >  r,  g(r)  >  0.  Thus 
there  exists  a  unique  solution  to  g(x)  ■  0,  say  6,  and  8  >  r.  But 
g(s)  <  0  iff 


-  >  iff  1  >  “I 

s  K(s)  n  ^  s+tj 


i 

r-»  *-*  elf 


n  ^  s+tj 


(2.5) 


i  ”  “~i 

therefore  g(s)  <  0.  Thus  the  unique  solution  8  is  such  that 


s  >  8  >  r. 

Now  that  we  know  the  unique  solution  of  g(x)  ■  0  exists  we  check 
that  it  is  indeed  the  maximum  likelihood  solution.  Since  by  substituting 
(2.3.1)  into  (2.2.3)  we  have 


1  3L 
n  38 


lr-Jl 


rs+8  -28r 


3  L 

it  is  sufficient  to  check  that  jjj- 


8-r 


K(8) 

”  °»  lei 


8- 


<  0.  By  substituting 


into  the  above  we  find  the  first  is  1/K(r)  >  0  and  the  second  is 
-1/s  +  1/K(s)  <  0  by  Equation  (4.5).  || 

We  now  present  two  methods  of  finding  8. 

Method  I:  If  the  data  satisfies  the  inequality 


2s  <  3r  +  min(t^t...,t  ) 


(2.6) 


then  the  Newton  iteration  procedure 


n+1 


B  S(B") 

n  '  8’<Sn) 


n-0,1, . . . 


will  converge  to  8  for  all  initial  points  r  <  8q  <  s ,  where  the 
function  g  was  defined  in  (2.1.3). 
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Proof.  To  assure  ourselves  that  the  Newton  method  will  converge 
It  Is  sufficient  that  g*,g"  do  not  vanish  for  r  <  x  <  s  since  they 
are  continuous. 

By  referring  to  the  formula  for  g'(x)  in  (2.4.1)  we  see  that 
for  r  <  x  <  s,  the  first  term  is  negative  since  K’(x)  >  1.  The 
second  term  is  x-r-K(x)  which  for  x  ■  r  is  negative  and  decreases 
for  x  >  r.  Thus  g'  (x)  <  0  for  r  <  x  <  s  and  cannot  vanish  in  that 
interval. 


We  now  show  that  g"(x)  >  0  for  r  <  x  <  s.  Note 


(2.6.1) 


g"(x)  -  2[1-K"(x)]  -  (x-r)K"(x), 


and  by  taking  derivatives  of  Equation  (2.4.2)  we  note 


K'.lxl  „  .1  }-2,2  _  — 1 —  .  }-3 

2K3(x)  n  U  nK(x)  U  V  * 

Thus  g"(x)  >  0  for  r  <  x  <  s  iff 


1-K*  (x)  <  (x-r)K11  (x) 
K3(x)  2K3(x) 


But  the  left-hand  side  of  (2.6.1)  Is  equal  to 


[n  2:(x+t1)"1][^E(x+ti)'2]. 


Multiplying  both  sides  of  (2.6.1)  by  (x-r)J  and  letting  X  (in  this 
argument  only)  denote  the  random  variable  which  takes  the  values 
(x-r)/(x+t^)  for  i>l,...,n  with  equal  probability,  we  find  that  (2.6.1) 
is  equivalent  with 


EX  EX2  -  (EX)3  >  EX  EX 3  -  (EX2)2. 


(2.6.2) 


Recalling  that  £n  EXV  is  convex  in  v  we  note  that  each  side  of 
the  inequality  (2.6.2)  is  positive. 

Since  EX  >  0,  we  can  divide  (2.6.2)  by  it  and  obtain 

EX2  -  (EX)2  >  EX3  -  . 

But  because  (EXV)^V  is  nondecreasing  in  v  we  know  that  (EX2)2  >_  (EX)^. 
Hence 

EX*  .  JSiii  <  ex3  .  (ex,  3 

and  it  is  sufficient  for  (2.6.2)  to  show  that  EX2  -  (EX)2  >_  EX'*  -  (EX)^, 

2  2 

or  equivalently  to  show  EX  (1-X)  2L  (EX)  (1-EX) .  But  for  any  convex  ip  we 

2  3 

know  EiKX)  >.  >KEX).  Thus  y(x)  ■  x  -  x  is  convex  for  x  <  1/3. 
Therefore,  a  sufficient  condition  that  gM(x)  not  vanish  on  r  <  x  <  s  is 

<  j  for  all  i-l,...,n.  (2.6.3) 

One  can  check  that  (2.6)  implies  (2.6.3)  in  the  interval  needed. || 

Method  II:  If  the  data  satisfies 

2r  >  s,  (2.7) 

then  for  H  ■  j  K,  where  K  was  defined  in  Equation  (2.1.2),  set 

A(x)  =  r  +  H(x)  -  /H2(x)-r(s-r)  (2. 7. 0.1) 

and  for  all  initial  points  r<BQ<s,  as  n  ”  the  iteration 
A(n)(80)  converges  to  g. 

Proof.  Solve  g(x)  *0  in  the  form 


x2  -  2x[r+H(x)]  +  r[s+2H(x)]  ■  0 


by  considering  H(x)  as  a  constant.  Then  using  the  quadratic  formula 
and  selecting  the  appropriate  root  we  find  the  solution  is  A(x)  as 
defined  in  (2. 7.0.1). 


In  order  to  assure  ourselves  that  the  values  A(x)  are  real  for 
x  >  r  we  need  to  show 


H(x)  >  /r(s-r)  .  (2.7.1) 

A  sufficient  condition  that  this  is  true  is  (2.7)  as  we  now  show. 

At  x  ■  r  (2.7.1)  is  equivalent  with 

1/4^7  >  *  l  -A-  , 

n  f  r+t^  * 

but  by  the  relationship  of  harmonic  and  arithmetic  means  we  know 


1(1  +  ±) 
2Kr  t ' 


and  hence 


Thus  by  assumption 
since  we  know  H  > 


i  Iy_L_<  I(i  + I)  .1  . 


H(r) 


r+tt  2  r  r' 


(2.7),  the  inequality  (2.7.1)  will  be  satisfied 
r  and  H'  >  0  for  all  x  >  0. 


One  notes  that 


A'  -  H' 


1  - 


H 


/H2-r(s-r) 


and  by  the  noted  properties  of  H  we  see  that  A'  <  0.  Hence  A  is 

monotone  decreasing.  Thus  if  we  can  show  for  r  <  <  s  we  have 

(2) 

A  (Bq)  bounded  (here  (2)  indicates  composition  with  itself)  then  by 
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Che  uniqueness  of  the  maximum  likelihood  estimate  (which  we  have  seen 

(n) 

previously)  and  the  completeness  of  the  real  line  we  have  Av  '(B^)  3. 

From  the  definition,  using  the  first  two  terms  of  the  binomial 
expansion  we  see, 

f  *  2H(Bpj  1  A(V  -  r  +  H(B0)>  (2.7.2) 

Applying  A  again  we  find 

A[r  +  fnoy  1  1a(2V  -A[r  +  H(80>1- 

But  by  (2.7.2)  we  have 

r  +  2H[r+H(BQ)]  -  a(  )(S0)  -  r  +  Hlr  +  2H(Bq)  (2.7.3) 

For  this  argument  let  Y  be  the  random  variable  taking  the  values 
t^  for  i«l,...,n  with  probability  1/n,  now  since  the  function  of 
y  defined  by  2/(x+y)  is  convex  in  y  and  EY  -  s  we  have  by 
Jensen's  inequality 

1  ,  rr  2  .  1  _  _L_ 

H(x)  *  ^x+Y;  -  x+EY  *  x+s  ' 

Applying  this  inequality  to  the  right-hand  side  of  (2.7.3): 

Using  (2.7.1)  we  obtain  a  condition  known  to  be  sufficient  for  convergence, 

r  <  A^(Bq)  <  s  +  2r  +  Vr(s-r)  .  |  | 

We  claim  that  the  conditions  (2.6)  and  (2.7)  are  not  stringent  in 
practice  and  will  nearly  always  be  satisfied. 
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3.  The  Mean  Mean 

Since  all  Iteration  methods  are  functionally  related  to  how 
close  the  initial  guess  is  to  the  answer  sought,  we  take  6q  “  6  where 

8  -  (SR)*5  (3.1) 

where  S  and  R,  now  considered  as  random  variables  and  written  in  the 
upper  case,  were  defined  in  (2.1.1).  We  choose  as  the  initial  estimate  for 
the  median  life  8  the  geometric  mean  of  the  harmonic  and  arithmetic  means 
of  the  sample  lives  in  the  data:  we  call  this  estimate  8  the  mean  mean. 

We  assert  that  the  great  utility  of  this  estimate  which  we  shall  demonstrate 
subsequently  can  be  accounted  for  in  the  following  comments. 

Theorem  3.1.  8  is  a  consistent  estimate  for  8* 

Proof.  By  (1.3)  and  the  strong  law  of  large  numbers,  with  probability 
one 

1  n  2 
sn  -  “  I  Tt  +  ET  -  8(1  +  ^-)  as  n  ->  ». 

11^1  !L 

But  also  ~  ■  —  £  (— )  and  —  ,  by  Theorem  1.1  has  the  distribution 
^  n  1  i  i 

F(a,  —)  so  that  with  probability  one 

P 

jp  -*■  ECj)  =  j(l  +  •Sj')  as  n  ->  <*>. 
n  p 

2  2 

Thus  clearly  (8  )  *  S^R^  ®  as  n  00  with  probability  one.|  | 

We  now  prove 

Theorem  3.2.  If  8  is  a  fixed  point  of  H  *  ^K,  then  8=6. 
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Proof.  Using  the  hypothesis  that  H(8)  *  8  we  find  from 
(2.1.3) 

g(B)  =  (8) 2  -  28r  -  2(B)2  +  (8)2  +  2r8  =  0 

but  by  definition  g(8)  *  0  implies  8=8,  since  8  is  the  unique 
solution  of  g(x)  ■  0  for  r  <  x  <  s.|  | 

A  V  A  - 

Theorem  3.3.  If  8  is  a  fixed  point  of  H  =  —  ,  then  8=8. 
Proof.  By  definition  8  is  a  fixed  point  of  A,  hence 

8  =  A(6)  =  r  +  H(8)  -  /h  (B)-r(s-r) 

but  by  hypothesis  H(8)  =8  so  that  8  =  r  +  8  -  / [8)2-r(s-r)  and 
solving  for  8  shows  that  8  =  ^rs.  |  | 

If  in  practice  the  conditions  should  obtain  for  which  the  mean 
mean  is  a  fixed  point  of  H,  within  a  reasonable  approximation,  it 
would  follow  from  the  preceding  two  results  that  within  that  degree 
of  approximation  we  could  take  the  mean  mean  as  the  maximum  likelihood 
estimate  of  8. 


We  now  present  some  sufficient  conditions  that  the  mean  mean  of 
a  set  of  numbers  be  a  fixed  point  of  H. 


Theorem  3.4.  If  t^,...^^  satisfy  the  relations 


t2i-l  =  3Ti’  fc2i  =  6y,Ti  for  i“l»***»k 


(3.2) 


or  the  relations 


t 


2i-l 


for  i=l , . . .  ,k 


(3.3) 


where  are  any  set  of  positive  numbers,  then  the  mean  mean 

8  from  these  samples  satisfies  exactly  the  identity  H ( B )  =  8. 


T-  * 
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Proof.  By  definition  we  must  show  that  1  *  2g/K(g)  but  note 
that  in  the  case  (3.2)  that  3=6  since 


s_ 

6 


_L_ 

2k 


Ti  + 


k 

I 


1 


_1_ 

rB 


K(B) 


1 

2 


Now  in  the  case  (3.3)  we  see  that  6*1  and  the  remainder  of  the 
argument  follows  from  that  just  given  with  g  *  1. 

Corollary  3.4.  The  conclusion  is  the  same  if 


in  case  (3.2)  we  add  t2k+1  =  8»  and 


in  case  (3.3)  we  add  t 


2k+l 


1. 


We  now  claim  that  any  sample  from  a  distribution  in  £7  tends 
to  act  like  the  sets  of  numbers  described  in  (3.2)  and  (3.3).  If 


T^,...,T2k  are  independent  random  variables  each  with  the  distribution 

F(a,3)  in  then  by  Theorem  1.1  we  see  that  , ...»  cp-  each  has 

1  k 

the  distribution  F(m,l)  which  is  exactly  the  same  as  the  distribution 


of 


k+1  2k 

-  ,...,  — ~  .  Thus  we  see  that  a  sample  tends  to  satisfy  (3.2). 


6 


6 


But  on  the  other  hand,  if  T.,,...,T,  ,  are  all 

1  k  k+1  X2k 

independent  and  identically  distributed  by  F(a,l),  then  8Tlt.. . ,gTk 


have  distribution  F(a,B)  while 


_L_ 

T  f  *  *  •  9  rp 

k+1  *2k 


also  have  the  same 


common  distribution.  Again  we  see  that  a  sample  tends  to  satisfy  (3.3). 


Thus  we  see  why  6  should  be  near  a  fixed  point  of  H. 


t 
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We  now  state 

Theorem  3.5.  If  T  has  the  distribution  F(a,l)  in  &,  then 
E(~)  -  1  for  all  a  >  0. 


Proof.  By  definition 


~  dF(t:a,l). 


Making  the  change  of  variable  y  *  1/t  in  the  third  integral  and 

realizing  that  F(—  :  ct,l)  *  1  -  F(s:a,l)  we  find  that  the  first  and 
s 

third  terms  cancel.  Hence 


ECjJf)  -  1  -  91(0)  -  \  .  1 1 

We  use  this  theorem  to  obtain  a  result  which  says  that  under  any 
conditions,  for  a  sufficiently  large  sample,  3  will  be  nearly  a  fixed 
point  of  H.  This  is  because  3/H(3)  approaches  unity  as  the  sample 
size  approaches  infinity.  Note  that  this  result  is  not  obtainable 
directly  from  Slutsky's  theorem,  given  e.g.,  p.  255,  Ref.  [4],  since  H 
itself  is  a  random  function  depending  upon  the  sample  size  despite  our 
suppression  of  this  fact  with  our  notation. 

Theorem  3.6.  The  random  variable  3_/H(3_)  converges  in  probability 
■  —  —  ■  n  n 

to  1  as  n  -*•  oo. 


Proof.  For  this  argument  let  us  set 


where  each  T  has  distribution  F(a,3).  Of  course  we  see  that 
S  (3  )  =  g  /H(3  ).  It  is  sufficient  to  show,  by  Theorem  3.5,  that 


Sn(6„)  5  S(B)  -  E(^jr). 


Now 
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P[|Sn(6n)-S<3)|  >  e]  i  P[|Sn(3n)-S(6n)|  >_c]  +  P[  |s(Bn)-S(B)  |  >.  e]. 


Since 


E(^?)2  *  f  - ^“7  dF(t:a,B/x)  1  1 

0  (1+t) 


1 

we  conclude  that  var(S  (x))  £  —  uniformly  in  x. 


By  Chebyschef f ’ s  inequality  as  n  -+  °° 


var  S  (B  )  C 

P[  |S  (3  )-S(8  )  |  >  e]  <  - 5—°-  <  -4-  o. 

nn  n  —  —  l  —  i. 

e  ne 


But,  of  course,  the  second  term  approaches  zero  by  the  consistency  of 

8n  for  8-11 

The  asymptotic  behavior  of  maximum  likelihood  estimates  such  as 
8  is  well  know.  We  now  prove  an  important  result  concerning  the 
asymptotic  distribution  of  the  mean  mean  8. 


Theorem  3.7.  For  n  sufficiently  large  the  estimator  3n  is 

asymptotically  distributed  by  F( —  ,8)  in  9  where 

Jn 


e2  =  (l  +  -^)/(l  +  4>2. 


(3.3.1) 


n  n  1 

Proof.  By  definition  8  *  where  Y  =  (  1  T  ) / (  ^  T  ) 

nn  n  1  j  1 

with  T^  independently  and  identically  distributed  by  F(a,l).  Let 


U.  =  Z./l  +  -y-  Z2 
li  4  l 


for 


i  =  l , . . . ,n 


(3.4) 


where  are  independent  s)?(0,l)  variates.  By  (1.5)  we  may  write 


Y  =  (1  +  ax  )/(l  -  ax  ) 
n  n  n 


(3.5) 


which  is  nearly  unity. 
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Ue  also  state 

Theorem  3.P.  The  estimate  6  is  for  each  fixed  n  biased 
- . . i  -  n 

by  a  factor  b^  which,  for  a  sufficiently  small,  is  given  by 

bn-i  +  f^+Ota4) 

with 

var(8  )  ■  +  O(a^) . 

n  n 

Proof.  The  proof  follows  using  similar  techniques  and  notation 
to  that  utilized  in  Theorem  3.7. 


For  n  fixed,  any  real  h  and  any  random  variables  i*l,...,n 

w*_  define  <U.  >  =  —  ]>  U .  .  The  bias  factor  of  8  is  b  *  e/T"",  where 
i  n  i  n  n  n’ 

Y  was  defined  in  (3.4).  From  (3.5)  we  have 
n 


-  <Ut>  I1  +  Y<V 


(3.8) 


where  here  the  U.  are  as  defined  in  (3.4).  On  the  set  B  =  H  [  I Z .  I  <  —  ], 
1  a  1  i 1  a 

we  may  apply  the  binomial  expansion  to  the  conditional  random  variable 
<U  >|B  =  l  (!s)(f)2j<Z^j+1>. 

1  °  jt0  J  2  1 

2  2 

Since  B^  C  [<Z^>  <  —  ]  by  using  (3.8)  we  have,  after  applying  the 

a 

binomial  expansion  and  simplifying, 


v  In  V  2k  v  72j  +  l  v2  M 
xnlB  =  2  a  >  c  <Z  J  ><Z  > 

k=0  j-0  J 

l,  _L  _  A  -!_[< 

where  for  notational  simplicity  we  set  c,  =  (_?)(,  2)4  J2J 

kJ  j  k-J 

definition  (3.5) 


From  the 


*/Y~  =  (1  +  nX  )'2(1  -  aX  )~''2 , 
n  n  n 
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again  since  c  M^nl  <  ”“1 »  by  expanding  and  simplifying  we  obtain 


By  rearranging  we  have 


where 


JT  |B  -  f  akc.Xk 
"  1  ~  k  n 


n  1  a 


k-0 


ck  *  S  (VhJ.k-d1- 

K  i=Q  1  K  1 


By  definition 

00 

b  »  S  akc,  E[X  IB  f]k  +  Z[JT  # BC|] 
n  k  n'  a*  1  n  1  a* 

where  fB|  is  the  indicator  function  of  that  set  B.  Note 


(3.9) 


«  <EYn)(l  -  P(Bb)] 

and  P(B  )  -  I'M—)  -  ?((—))"  *  1  as  a  -»  0. 
a  a  a 

Using  (3.7) 

X  |B  -  <Z,>  +  cx2cin<zJ><Z,>  + a2c11<Z.3>  +  0(a4) 
n1  a  i  10  i  i  Hi 

[XjB  J2  -  <Zt>2  +  2a2c10<Z2xZi>2  +  2a2c11<Z3xZi>  +  0(a4) 
[XjBj3  =  <V3  +  0(a2) 
and  for  k  >_  4 

[XjB/  -  0(ak"1). 
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/l 

Taking  the  expectation  and  realizing  that  | [ | |  <  —  ]  is  a 
symmetric  random  variable  we  see 

EIXJBJ]  -  OCa4) 

EIXJ\»)2  -  ^+0(“2> 

EIX„?BJ13  '  0(“2)’ 

Since  c«  ■  ■  1,  c2  “  c3  “  2^  we  conclu<^e  from  (3.8)  that  (3. A) 

holds. 

Consider  the  variance  of  8  .  We  know 

n 

var(0  )  *=  g2(EY  -b2). 
n  n  n 

Proceeding  as  before  have 

Y  I B  *  1  +  2aX  +  2a2X2  +  2a3X3  +••• 
n1  a  n  n  n 

EY  tB  I  -  1  +  ~  +  0(a4) 
n*  a'  n 

and  by  squaring  (3. A)  and  subtracting  we  obtain  (3.5). || 

This  completes  our  theoretical  study  of  the  behavior  of  the  mean 

mean.  In  summary,  we  have  seen  that  8n  is  a  consistent  estimator  of 

B  which  for  a  small,  say  less  than  1/2,  is  almost  unbiased  while 

2 

having  a  variance  of  approximately  (aB)  /n.  Moreover,  we  claim  that 
under  this  condition,  which  we  shall  later  empirically  verify,  Br  is 
virtually  the  maximum  likelihood  estimator  whose  optimal  properties 


are  well  known. 


I 
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A.  Numerical  Examples  and  Comparisons 

Let  us  first  perform  an  empirical  comparison  of  the  two  iterative 
procedures  for  obtaining  the  maximum  likelihood  estimator  g,  and  check 
to  see  if  the  known  consistency  of  this  estimate  is  of  practical 
significance.  We  shall  obtain  a  set  of  observations  from  the  random 
variable  defined  by  (1.5)  by  selecting  a, 3  and  then  using  a  library 
machine  program  for  the  IBM  360  to  generate  the  pseudo-random  normal 
observations.  We  know  that  the  computed  random  variable  has  a  distribution 
in  0. 

For  given  t^,...,t  generated  by  this  method  we  compute  the  sample 
mean  mean  6  *  /sr,  where  s  and  r  are  defined  in  (2.1.1).  (As  we  know 
3  is  itself  a  consistent  estimate  of  6.)  We  stop  the  iteration  when  the 
successive  iterates  agree  within  a  prescribed  e. 

The  results  of  this  simulation  are  presented  in  Table  I.  Method  II, 
which  computationally  is  the  simpler,  appears  to  work  as  well  for  a  <_  .5 
as  the  well-known  Newton  procedure,  however,  it  does  not  work  at  all  for 
values  of  a  as  large  as  2. 

The  important  point  is  that  for  a  <  .5  there  is  no  real  need  to 
compute  anything  but  3  since  both  estimates  agreed  to  within  three 
significant  digits.  It  is  doubtful  that  any  further  iteration  by  either 
method  would  really  improve  the  accuracy  of  the  estimate. 

The  computational  simplicity  of  the  estimate  S,  as  well  as  the 
advantages  of  its  asymptotic  distribution  being  in  the  class  not  to 

mention  the  properties  which  were  set  out  in  the  last  paragraph  of  the 
proceeding  section,  can  all  be  utilized  if  the  range  of  a  which  is 


V 
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encountered  in  practice  is  no  larger  than  say  1/2.  Thus  in  fatigue 
studies  it  would  be  desirable  that  the  appropriate  value  of  a  were 
within  this  range. 


I 
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Table  I 

Comparison  of  methods  to  find  MLE  within  e  for  g  =  100,  e  =  .001 

sai?ple  number  of  iterations 

slze  _  for  convergence  within  e 


n 

6 

6 

a 

I 

II 

5 

105.119 

105.117 

.184 

2 

2 

10 

99.611 

99.580 

.358 

2 

3 

25 

87.279 

87.281 

.602 

2 

2 

5  50 

110.117 

110.125 

.497 

2 

2 

100 

101.201 

101.194 

.590 

2 

2 

250 

99.676 

99.702 

.496 

2 

3 

500 

100.215 

100.207 

.476 

2 

2 

5 

79.180 

79.200 

.363 

2 

2 

10 

67.302 

68.584 

.986 

3 

8 

25 

113.710 

116.176 

1.095 

3 

10 

50 

86.470 

86.225 

1.091 

2 

7 

100 

113.666 

114.202 

.866 

2 

6 

250 

105.-968 

106.156 

.964 

2 

6 

_500 

94.973 

95.042 

1.016 

2 

5 

5 

305.983 

310.994 

1.748 

3 

* 

10 

48.948 

52.314 

1.654 

3 

* 

25 

125.3'’? 

131.245 

1.796 

3 

* 

50 

96.638 

99.622 

1.652 

3 

* 

100 

87.675 

88.032 

2.110 

2 

* 

250 

99.848 

99.634 

2.030 

2 

* 

_ 500 

96.640 

96.696 

2.127 

2 

* 

*Method  II  failed  to  be  applicable  since  for  some  point  6  in  the 
iteration  we  have  K(6n)  <  r(s-r)  and  6n+1  is  no  longer  real. 


-23- 


We  now  confront  this  family  £7  of  distributions  with  some  actual 
fatigue  data.  We  choose  some  extensive  data  on  the  fatigue  life  of 
6061-T6  aluminum  coupons  cut  parallel  to  the  direction  of  rolling  and 
oscillated  at  18  cycles  per  second.  Some  of  this  data  has  been  reported 
earlier  in  [2], 


For  a  maximum  stress  per  cycle  of  31,000  psi  we  give  the  101 

-3 

observations  of  lifetimes  in  cycles  X10 

SAMPLE  SIZE  =  101 


70 

90 

96 

97 

99 

100 

103 

104 

104 

105 

107 

108 

108 

108 

109 

109 

112 

112 

113 

114 

114 

114 

116 

119 

120 

120 

120 

121 

121 

123 

124 

124 

124 

124 

124 

128 

128 

129 

129 

130 

130 

130 

131 

131 

131 

131 

131 

132 

132 

132 

133 

134 

134 

134 

134 

134 

136 

136 

137 

138 

138 

138 

139 

139 

141 

141 

142 

142 

142 

142 

142 

142 

144 

144 

145 

146 

148 

148 

149 

151 

151 

152 

155 

156 

157 

157 

157 

157 

158 

159 

162 

163 

163 

164 

166 

166 

168 

170 

174 

196 

212 

Setting  t  =  .0001  we  find  using  Method  I,  starting  with  fj  =  131.819454 
that  in  three  iterations 


6  =  131.81903 
\  =  .17037302, 


while  with  Method  II,  with  the  same  value  of  8  as  initial  guess,  we 


find  that  in  two  iterations 


6  =  131.81895 
a  =  .17037022. 


« 


From  the  graphical  comparison  of  the  fitted  dis  ribution  to  the 

empiric  cumulative  in  Figure  1,  we  see  that  9  provides  an  adequate 
explanation. 


Cycles  X  10'3 


Figure  1 

The  empiric  cumulative  and  the  distribution  F(a,B)  for  fatigue  life 
at  a  stress  of  31,000  psi. 


*  • 
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For  a  maximum  stress  of  26,000  psi  we  give  the  102  observations 

-3 

of  lifetimes  in  cycles  X10 
SAMPLE  SIZE  =  102 


233 

258 

268 

276 

290 

310 

312 

315 

318 

321 

321 

329 

335 

336 

338 

338 

342 

342 

342 

344 

349 

350 

350 

351 

351 

352 

352 

356 

358 

358 

360 

362 

363 

366 

367 

370 

370 

372 

372 

374 

375 

376 

379 

379 

380 

382 

389 

389 

395 

396 

400 

400 

400 

403 

404 

406 

408 

408 

410 

412 

414 

416 

416 

416 

420 

422 

423 

426 

428 

432 

432 

433 

433 

437 

438 

439 

439 

443 

445 

445 

452 

456 

456 

460 

464 

466 

468 

470 

470 

473 

474 

476 

476 

486 

488 

489 

490 

491 

503 

517 

540 

560 

Again  choosing  e  *  .0001  we  find  using  Method  I  starting  with 
6  *  392.765189,  that  in  two  iterations 

B  =  392.76367 
a  -  . 16141957 

and  that  for  Method  II  in  two  iterations 

6  =  392.76416 
a  =  .1614195. 


See  Figure  2  for  a  graphical  comparison 


For  a  maximum  stress  of  21,000  psi  we  have  101  observations  of 
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lifetimes  in  cycles  X10 


SAMPLE  SIZE  =  101 


370 

706 

716 

746 

785 

797 

844 

855 

858 

886 

886 

930 

960 

988 

990 

1000 

1010 

1016 

1018 

1020 

1055 

1085 

1102 

1102 

1108 

1115 

1120 

1134 

1140 

1199 

1200 

1200 

1203 

1222 

1235 

1238 

1252 

1258 

1262 

1269 

1270 

1290 

1293 

1300 

1310 

1313 

1315 

1330 

1355 

1390 

1416 

1419 

1420 

1420 

1450 

1452 

1475 

1478 

1481 

1485 

1502 

1505 

1513 

1522 

1522 

1530 

1540 

1560 

1567 

1578 

1594 

1602 

1604 

1608 

1630 

1642 

1674 

1730 

1750 

1750 

1763 

1768 

1781 

1782 

1792 

1820 

1868 

1881 

1890 

1893 

1895 

1910 

1923 

1940 

1945 

2023 

2100 

2130 

2215 

2268 

2440 

With 

e  *  .0001, 

8  -  1336 

.56547 

we  find 

that  in 

fifty-one 

Iterations 

Method  I  failed  to  converge  yielding  as  values  at  that  time 


8  =  1336.3779 
a  -  .31029648, 


but  that  Method  II  converged  in  three  iterations  to 


8  -  1336.3784 
a  *  . 31029648. 


We  believe  that  the  extreme  stringency  of  e  and  round-off  error 
in  the  machine  arithmetic  caused  Method  I  to  fail  to  converge  and  not  a 
theoretical  deficiency  in  the  method  itself. 

We  can  tentatively  conclude  from  this  evidence  that  in  fatigue 
applications  the  appropriate  range  of  a  is  sufficiently  small  as  to 
allow  the  use  of  8  as  an  estimate  of  the  median  8  and  thus  one  can 


utilize  the  properties  which  we  have  previously  discussed. 
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